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Abstract 

Vibrational energy transfer of the amide I mode of N-methylacetamide (NMA) 
is studied theoretically using the vibrational configuration interaction method. A 
quartic force field of NMA is constructed at the B3LYP/6-31G+(d) level of theory 
and its accuarcy is checked by comparing the resulting anharmonic frequencies with 
available theoretical and experimental values. Quantum dynamics calculations for 
the amide I mode excitation clarify the dominant energy transfer pathways, which 
sensitively depend on the anharmonic couplings among vibrational modes. A ratio 
of the anharmonic coupling to the frequency mismatch is employed to predict and 
interpret the dominant energy flow pathways. 
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1 Introduction 



The dynamics of proteins involves high frequency vibrational modes that be- 
have quantum mechanically even at room temperature (T ~ 300 K corre- 
sponds to ~ 200 cm~^). One such class of vibrations is the amide I mode, 
centered around 1650 cm~^ and prominent in IR experiments due to its large 
oscillator strength. The amide I modes have been intensively studied by 2D-IR 
spectroscopy, which has provided new insights into their anharmonic couplings 
to the other vibrational modes and their intrinsic mode anharmonicities [1,2,3]. 
Theoretical studies have also been stimulated to characterize the details of the 
vibrational states [4,5,6,7], and to probe the origin of the inhomogeneous line- 
shape broadening of the amide I mode [5,8,9,10]. 



In this study, we focus on the vibrational energy relaxation (VER) process of 
the amide I mode. Recent pump-probe experiments [11,12] have found that 
the VER of the amide I mode occurs on a a sub-picosecond timescale after 
its excitation. Nguyen and Stock [13] have studied this phenomenon using 
molecular dynamics and instantaneous normal mode analysis. Their results 
were found to be in qualitative agreement with the experimental results. Nev- 
ertheless a more precise quantum mechanical description of the VER process 
is desirable. Two of the authors [14] have recently proposed a time-dependent 
perturbation theory to describe VER, and applied the formula to a small 
peptide- like molecule, N-methylacetamide (NMA), in heavy water using an 
empirical force field. They also observed VER on a time scale (~ 0.5 ps) sim- 
ilar to that obtained experimentally, and further proposed a mechanism for 
the VER process. 



The present study explores the VER processes of NMA more accurately by 
directly solving the Schrodinger equation for molecular vibrations on ah ini- 
tio potential energy surface (PES). The NMA molecule in vacuum is exam- 
ined, excluding water molecules included in the previous studies [13,14]. Quan- 
tum dynamics calculations are carried out by the vibrational configuration- 
interaction (VCI) method [15,16,17,18]. 24 vibrational modes (out of 30) are 
explicitly treated in the dynamics, excluding the 6 lowest-lying modes includ- 
ing the rotational motions of the two methyl groups. We note that while the 
multiconfiguration time-dependent Hartree (MCTDH) method [19] provides 
a more fiexible framework the VCI method is sufficient for the present aim to 
investigate the vibrational motion of a semi-rigid molecule. 
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2 Methods 



2. 1 Quartic force field 



The PES is approximated using a fourth-order Taylor series expansion around 
the equihbrium geometry, which is called quartic force field (QFF), as 

^{{Qi}) — -^'^^^kQk + ^ X] ^klmQkQlQm + ^ XI '^klmnQ kQ iQ mQ m (l) 
k ' k,l,m ' k,l,m,n 

where Qk and Uk denote the A:th normal coordinate and the associated har- 
monic frequency. The coefficients, tkim and Ukimn, are the third- and fourth- 
order derivatives of the PES. The above QFF can be recast in the form of the 
n-mode couphng representation (nMR) [20,21] as 



V{{Qi}) ~ Vi{Q^}) = y^^^ + y^Mil + ySMR^ ^2) 

V'^'^'m}) = i: (^ujlQl + y,,,Ql + ^u,,,,Qt^ , (3) 

v'^'^'aQ,}) = Y: {^^hkiQlQi + \ukkiiQlQ] + \y^kkkiQlQ^ , (4) 

yZMR(^^Q^^ = X! [iTy^klmQkQlQm + -T:^kklmQ\QlQm \ ■ (5) 
k,l,m ^ 

It is known that 3MR-QFF is enough to characterize the anharmonicity of 
molecules studied [20,21], we thus have neglected the fourth-order terms in- 
cluding Ukimn with k ^ I ^ m ^ n. 

In addition to the above 3MR-QFF, we study two types of approximate QFF: 
(1) 2MR-QFF that neglects V^^^ [4] and (2) partial-QFF that takes into 
account only the terms associated with the amide I mode, that is, tkim — 
and Ukkim = if the subscripts do not include the amide I mode. Note 
that the previous perturbation calculation [14] also employed partial-QFF to 
calculate the reduced density matrix. By adding the normal mode kinetic 
energy K — Z^fe-Pfe/2, we have the full approximate vibrational Hamiltonian 
H + V{{Qi}) for the system. 



2.2 Vibrational CI method 



The vibrational self-consistent field (VSCF) calculation is first carried out 
for the vibrational ground state to obtain the basis functions for the VCI 
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calculations. The VSCF wavefunction is expressed as a direct product of one- 
mode functions or modals as 



K'^'^-h^m. (6) 
1=1 

where n and / denote the vibrational quantum numbers and the degrees of 
freedom, respectively. The modals are determined by 



2QQ2 >xx.-.. 



(7) 



This VSCF equation is solved for the vibrational ground state (n=0), and 
the virtual modals constitute the VSCF configurations to be used in the VCI 
calculations. The VCI wavefunction is expressed as a linear combination of 
VSCF configurations as 



*r = ECmn*I'^^. (8) 
m 

The VCI wavefunction and energy levels are obtained by diagonalization of 
the VCI matrix 

H^^^i^-'j^^m^i^^^). (9) 



In this study, the modals were expanded in terms of the harmonic oscillator 
(HO) wavefunction. The number of HO wavefunctions employed were 11, 9, 
7, and 5 for {0(7), {^O). 0(i2)}^ |^(i3)_^(23)|^ {0(24).0(3O)}^ respec- 

tively. The mode index is labeled in the increasing order of the frequency. See 
Table 1. The 6 lowest-lying modes were kept frozen. The VSCF configurations 
were constructed by allowing the excitation up to 10 quantum numbers, and 
selected in the increasing order of the energy until the upper limit of the VCI 
space, denoted Nqi^ was achieved. The VSCF/ VCI calculations were carried 
out using the SINDO code [22] for non-rotating molecules. 



2.3 Quantum dynamics 



Once the eigenvalues {{E-a}) and the eigenfunctions ({^'n*"^}) are obtained, it 
is straighforward to carry out an approximate quantum dynamics simulation 
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Table 1 

Calculated harmonic (H.O) and anharmonic (VCI and PT2) frequencies of NMA 
based on the B3LYP/6-31G+(d) level of theory. The previous cc-VSCF and VCI 
results and the experimental results are also hsted for comparison. Units in cm~^. 



Mode 


H.0.« 


VCI^ 


PT2^ 


cc-VSCF^ 


VCI*^ 


Exp.^ 


7 


623 


619 


614 


636 


633 


619 


8 


630 


625 


613 


637 


685 


658 


9 


879 


869 


861 


891 


886 


857 


10 


1001 


983 


982 


1022 


1005 


980 


11 


1066 


1046 


1042 


1083 


1061 


1037 


12 


1107 


1093 


1080 


1119 


1099 


1089 


13 


1163 


1141 


1126 


1184 


1167 


— 


14 


1199 


1186 


1146 


1214 


1199 


1168 


15 


1292 


1272 


1256 


1283 


1253 


1266 


16 


1421 


1397 


1387 


1421 


1388 


1370 


17 


1473 


1448 


1432 


1468 


1442 


1419 


18 


1495 


1451 


1481 


1515 


1467 


1432 


19 


1500 


1453 


1459 


1557 


1483 


1432 


20 


1507 


1469 


1476 


1541 


1481 


1446 


21 


1529 


1495 


1481 


1566 


1505 


1472 


22 


1560 


1537 


1505 


1547 


1519 


1511 


23 


1751 


1725 


1725 


1751 


1727 


1707 


24 


3059 


3017 


2906 


2939 


2940 


2915 


25 


3060 


2996 


2940 


2940 


2983 


2958 


26 


3117 


3077 


2963 


2979 


2995 


2973 


27 


3132 


3093 


2977 


3014 


3043 


3008 


28 


3136 


3100 


2990 


2985 


3025 


2973 


29 


3149 


3077 


2990 


2993 


3056 


3008 


30 


3643 


3479 


3466 


3523 


3544 


3498 



" Harmonic frequencies at the B3LYP/6-31+G(d) level. 

* Based on 3MR-PES at the B3LYP/6-31+G(d) level. 

" Reference [4]. Based on 2MR-PES at the MP2/DZP level. 

Reference [7]. Based on partial-3MR-PES at the MP2/aug-cc-pVTZ level. 
^ Reference [24]. 
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n 

If we take the initial state to be the VSCF configuration, \E'(0) = ^X^*^^, then 
the overlap integral between ^{t) and ^y^*-^^ is calculated as 

= (*7SCF|^(^)^ ^ ^ CjnCnie-^^"*/'^. (11) 
n 

The absolute square of Oj{t) gives the probability of j at time t, Pj(t) — 
|Oj(t)p. Dynamics calculations were carried out with the initial VSCF con- 
figurations corresponding to the fundamental and first overtone of the amide 
I mode. We note that we also examined different initial states which are the 
superposition of states near the amide I mode fundamental or overtone, but 
the results do not severely depend on how to prepare the initial states [23]. 



3 Results 

3.1 Anharmonic frequency calculations of NMA: Accuracy of the PES 

Here we shall examine the accuracy of the PES derived from B3LYP/6- 
31-|-G(d) employing anharmonic frequency calculations because the vibra- 
tional quantum dynamics sensitively depends on anharmonicity of a system. 
We compare our result with the previous theoretical values obtained by cc- 
VSCF method [4] and very recently by VCl method [7] using MULTIMODE 
[17] as well as with the experimental values [24]. We have also computed 
the fundamental frequencies of NMA by the second-order perturbation theory 
(PT2) [25,26] using Gaussian03 [27]. Table 1 shows the resulting harmonic 
(H.O.) and anharmonic frequencies (VCl and PT2) based on the B3LYP/6- 
31+G(d) level of theory, together with the previous cc-VSCF results based on 
the MP2/DZP [4] and the VCl results based on the MP2/aug-cc-pVTZ [7]. 
The experimental results are taken from infrared spectra of NMA in nitrogen 
matrix [24]. Although the experiment was not in the gas phase, we expect 
that the nitrogen matrix should affect minimally for vibrational frequencies, 
and we regard these correct values. As is well known, the harmonic frequencies 
overestimate the correct frequencies but the anharmonic frequencies are closer 
to the latter. A close inspection shows that the PT2 result is the most accurate 
but marginally and the other three methods are almost comparable except for 
six CH stretching modes (~ 3000 cm^^), where our result seems to be crudest. 
However, we are interested in the amide 1 mode, whose frequency is around 
1700 cm~^, and the VER pathways from the amide I mode are usually toward 
lower frequency modes, in this paper we accept this level of accuracy. 
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3.2 Vibrational quantum dynamics of NMA-D: Fermi resonance 



To compare with the previous dynamics calculation using a force field [14] and 
time-resolved experiments [11,12], we deuterated NMA into NMA-D (CH3- 
ND-CO-CH3). We show the VER dynamics of NMA-D calculated by VCI 
method in Fig. 1. We consider two initial states: v = 1 (fundamental) and v = 
2 (overtone) excitation of the amide I mode in the VSCF base. In both cases, 
the initial "decay" appears to occur on a sub picosecond timescale, which is in 
accord with the previous studies [12,14]. However, the long-time dynamics (not 
shown) are quasi-periodic with resonant energy "transfer" between the amide I 
mode and the other modes. In order to induce irreversible decay, it is essential 
to include more bath modes in the form of solvent molecules so that energy 
can dissipate [14] . A vibrationally excited molecule with sufficient quanta can 
decay irreversibly even in a vacuum [28]. The present result indicates that 
the intramolecular energy "transfer" process is responsible for the initial fast 
"decay" (< 1 ps). 





Fig. 1. Quantum population dynamics of the vibrational excited states, v = 1 
(left) and v = 2 (right) of the amide I mode, on 3MR-QFF derived from 
B3LYP/6-31G+(d) method. Nci is set to 6000. The dominant VER pathways are 
also summarized in Table 2. 

The dominant VER pathways are described in Table 2, and the resonant bath 
modes are depicted in Fig. 2. We interpret these dominant pathways using the 
following Fermi resonance parameter [29] 



mv\f) 

AE 



(12) 



where \i) and |/) are the initial and final harmonic states, AV = V — 
J^k^lQlf^ anharmonic potential energy, and AE is the energy dif- 

ference between \i) and |/). For example, rj for a transition between \i) = \Si) 
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Table 2 

The dominant VER pathways in NMA-D and the corresponding Fermi resonance 
parameters when we excite the fundamental of the amide I mode. \mn) denotes 
the nth excited states of the mth mode, and all other modes are on the vibrational 
ground state. The 23rd mode is the amide I mode in this paper. The result is derived 
from the B3LYP/6-31G+(d)/3MR-PES. 



VER pathway Fermi resonance parameter 



|23i) ^ |92> 
|23i) ^ |12i7i) 

|23i) ^ IIO2) 
|23i) ^ |10i9i) 




0.082 
0.024 
0.014 
0.013 





Fig. 2. Normal mode vectors of the resonant bath modes for the amide I mode 
in deuterated NMA. (a) 7th (619 cm"^), (b) 9th (869 cm^^), and (c) 12th (1066 

cm" 



' mode. 



(the system mode is singly excited) and |/) = \kili) (two bath modes are 
singly excited) is evaluated as 



{Si\2tskkQsQkQi\kih) 



6h{uJs -ujk- uji) 



^Skl 



'A. A. A.(i3) 



Note that the Fermi resonance parameter was a key ingredient of the time- 
dependent perturbation theory in describing the reduced density matrix for 
the amide I mode [14]. Similar analysis has been applied to protein dynamics 
of myoglobin in terms of classical mechanics [30]. 

Table 2 also shows that the dominant VER pathways are well characterized 
by the above parameters for the v = 1 excitation. Figure 3 shows that other 
minor pathways take on values of r] that are smaller than ~ 0.01 and that 
the frequency matching condition is important to predict the Fermi resonance 
parameter. It is interesting to note that the dominant pathways (Fig. 2) are 
similar to those suggested in [11]. In the case oi v = 2 excitation (Fig. 1, 
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right), it is difficult to assign such a parameter because the relaxation processes 
can accompany higher order processes, but the dominant pathways such as 
1282) — >■ |23i, 12i,9i) and |232) —>■ \9i) are the frequency matching ones. We 
conclude that the Fermi resonance parameter or frequency matching condition 
is very important to characterize the VER pathways, and its consequence is 
discussed below. 

3.3 Discussions 

From these observations, we conchide that a relatively small number of reso- 
nant bath modes plays an essential role in the quantum dynamics. Simplifica- 
tion of the PES or a reduction in the number of modes, including the resonant 
bath modes, can result in an inaccurate estimate of the dynamics. This is 
illustrated in Fig. 4: the 3MR-QFF and partial-QFF results agree rather well 
up to 0.2 (0.5) ps for f = 1 {v = 2), while 2MR-QFF results deviate from the 
other two at the initial stage. This is because the 2MR representation misses 
resonances mediated through the three mode interactions. It is interesting that 
the vibrational frequencies can be calculated accurately using 2MR-QFF; the 
fundamental (first overtone) of the amide I mode is obtained as 1713 (3414), 
and 1716 (3421) cm^^ by the 2MR- and 3MR-QFF, respectively (the differ- 
ence is 0.2 %). However, the VER processes sensitively depend on the 3MR 
terms. 

We showed that our PES is accurate enough (Sec. 3.1), which is in accord 
with experiment [24] and comparable to other theoretical methods [27,4,7]. 
However, we note that the VER dynamics are much more sensitive to accuracy 
than the anharmonic frequency calculations. For example, Gerber's group did 
not report any strong resonant interaction between the amide I mode and 
other modes using MP2/DZP level of theory [4]. This is likely due to the level 
of theory for the PES (MP2/DZP), but another possibihty is due to the 2MR- 
PES in their calculations. As we showed above, 3MR-PES and 2MR-PES can 
provide similar anharmonic frequencies, but some resonance conditions should 
be missed in the latter, affecting the dynamics calculations. Importantly, 2D- 
IR spectroscopy can directly detect anharmonic coupling in a molecule [1,2,3], 
and such data may be utilized to select an appropriate level of methods for a 
particular molecule. 



4 Concluding remcirks 

Using the VCI method, we investigated the quantum dynamics of deuter- 
ated N-methylacetamide in vacuum. We demonstrated the applicability of the 
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1/|(0s -cOj -cOjl (cm) 




Fig. 3. Top: The Fermi resonance parameter as a function of the mode index (i, j). 
Middle: The 3rd order coupling strength \tsij \ as a function of the mode index (i, j). 
Bottom: The inverse of the frequency mismatch l/\iJs — — Wj| as a function of 
the mode index (i,j). 

method and were able to gain insight into the fundamental nature of VER 
in the molecule, relevant to the interpretation of IR and 2D-IR spectroscopy 
used as a probe of protein dynamics. The accuracy of the PES employed 
(B3LYP/6-31+G(d)) was checked by the anharmonic frequency calculations 
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0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 



t (ps) t (ps) 

Fig. 4. Quantum population dynamics of the vibrational excited states, v = 1 (left) 
and f = 2 (right) of the amide I mode, on 3MR-QFF, 2MR-QFF, and partial-QFF, 
derived from B3LYP/6-31G+(d) method. Nci is set to 6000. 

and by comparing with experiment and other theoretical methods. We clari- 
fied the energy fiow pathways from the v = 1 and 2 excitations of the amide 
I mode, and interpreted our results using Fermi resonance parameters, which 
represent the effective coupling strength betweem vibrational modes. This ap- 
proach will be extended to condensed phase systems by invoking QM/MM 
methods [31,32] or multiresolution methods [6,33,34] to deal with dephasing 
problems and to simulate 2D-IR signals [35]. 
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